Assessments of Vascular Permeability For Biomedical Imaging Studies

ABSTRACT

Described herein are methods and systems for analyzing biomedical images using new models. Example models include a linear reference region model and a reference agent model. In one example aspect, a computer-implemented method is provided. The method may involve determining, based on a set of biomedical images, a first concentration-activity curve and a second concentration activity-curve. Additionally, the method may further include determining a value of at least one pharmacokinetic (PK) parameter based on the first concentration-activity curve and the second concentration-activity curve and a linear model that relates the first concentration-activity curve to the second concentration-activity curve. The value of the at least one PK parameter may be determined based on application of a linear least square fitting algorithm to the linear model. Also, the method may include causing a graphical display to provide a visual indication of the value of the at least one PK parameter.

RELATED APPLICATIONS

This application claims priority to U.S. Provisional Patent ApplicationSer. No. 61/633,549 filed Feb. 13, 2012, entitled Improved Assessmentsof Vascular Permeability Using a Linear Reference Region Model forDynamic Contrast Enhancement Imaging Studies and claims priority to U.S.Provisional Patent Application Ser. No. 61/634,589 filed Mar. 2, 2012,entitled Improved Assessments of Vascular Permeability Using a ReferenceAgent Model for Dynamic Contrast Enhancement Imaging Studies, each ofwhich is herein incorporated by reference in its entirety.

STATEMENT OF U.S. GOVERNMENT INTEREST

This invention was made with government support under P30 CA023074awarded by the National Cancer Institute at the National Institutes ofHealth, and W81XWH-10-1-0188 awarded by the U.S. Army Medical Researchand Material Command. The government has certain rights in theinvention.

FIELD OF THE INVENTION

The disclosure herein relates generally to processing biomedical imagesand, in particular, to methods and systems for assessing vascularpermeability and flow using biomedical images.

BACKGROUND

Biomedical imaging techniques can be used to assess physiologicalparameters of a tumor, organ, or other tissue of interest. For instance,biomedical imaging techniques can be used to estimate vascularpermeability, the movement of fluids and molecules between vascular andextravascular compartments. Vascular permeability may characterize theability of a blood vessel wall to allow for the flow of molecules (e.g.,ions, nutrients, or water) or cells into or out of the blood vessel.Other examples of uses for biomedical imaging techniques may exist.

Current biomedical imaging techniques for assessing vascularpermeability and flow are inadequate.

SUMMARY

Described herein are new methods and systems for analyzing biomedicalimages. In one example, methods and systems that may provide forimproved speed of analysis and less sensitivity to both noise andtemporal resolution are described. Such an example may enable theanalysis of biomedical images obtained in a clinical setting. In anotherexample, methods and systems are provided that involve estimation of therelative permeability and/or flow of two contrast agents within the sametissue of interest, advantageously eliminating some of the physiologicalvariables that might otherwise affect the analysis.

In one example aspect, a computer-implemented method is provided. Themethod may involve determining, based on a set of biomedical images, afirst concentration-activity curve and a second concentrationactivity-curve. Each of the first concentration-activity curve and thesecond concentration-activity curve may indicate a concentration of arespective contrast agent within a respective region of tissue.Additionally, the method may further include determining a value of atleast one pharmacokinetic (PK) parameter based on the firstconcentration-activity curve and the second concentration-activitycurve. Determining the value of the at least one PK parameter mayinclude (i) determining a linear model, including the at least one PKparameter, that relates the first concentration-activity curve to thesecond concentration-activity curve and (ii) determining the value ofthe at least one PK parameter based on application of a linear leastsquare fitting (LLSQ) algorithm to the linear model. Also, the methodmay include causing a graphical display to provide a visual indicationof the value of the at least one PK parameter.

In another example aspect, a non-transitory computer-readable mediumhaving instructions stored thereon is provided. The instructions mayinclude determining, based on a set of biomedical images, a firstconcentration-activity curve and a second concentration activity-curve.Each of the first concentration-activity curve and the secondconcentration-activity curve may indicate a concentration of arespective contrast agent within a respective region of tissue.Additionally, the instructions may further include determining a valueof at least one pharmacokinetic (PK) parameter based on the firstconcentration-activity curve and the second concentration-activitycurve. Determining the value of the at least one PK parameter mayinclude (i) determining a linear model, including the at least one PKparameter, that relates the first concentration-activity curve to thesecond concentration-activity curve and (ii) determining the value ofthe at least one PK parameter based on application of a linear leastsquare fitting (LLSQ) algorithm to the linear model. Also, theinstructions may include causing a graphical display to provide a visualindication of the value of the at least one PK parameter.

In a further aspect, a system comprising at least one processor, anon-transitory computer-readable medium, and program instructions storedon the non-transitory computer-readable medium is provided. Theinstructions may be executable by the at least one processor todetermine, based on a set of biomedical images, a firstconcentration-activity curve and a second concentration-activity curve.Each of the first concentration-activity curve and the secondconcentration-activity curve may indicate a concentration of arespective contrast agent within a respective region of tissue.Additionally, the instructions may be executable to determine a value ofat least one pharmacokinetic (PK) parameter based on the firstconcentration-activity curve and the second concentration-activitycurve. Determining the value of the at least one PK parameter mayinclude (i) determining a linear model, including the at least one PKparameter, that relates the first concentration-activity curve to thesecond concentration-activity curve and (ii) determining the value ofthe at least one PK parameter based on application of a linear leastsquare fitting (LLSQ) algorithm to the linear model.

All embodiments of the invention disclosed herein can be combined withother embodiments or combinations of embodiments unless the contextclearly dictates otherwise.

These as well as other aspects, advantages, and alternatives, willbecome apparent to those of ordinary skill in the art by reading thefollowing detailed description, with reference where appropriate to theaccompanying drawings.

BRIEF DESCRIPTION OF THE FIGURES

FIG. 1 is a flow chart of an example method for determining a value ofone or more pharmacokinetic parameters.

FIG. 2 is a conceptual illustration of an example model used for thequantitative analysis of biomedical image data.

FIG. 3 is a conceptual illustration of another example model used forthe quantitative analysis of biomedical image data.

FIGS. 4-7 show the results of example computer simulations performed toevaluate the use of a linear model.

FIGS. 8-9 show the results of an example experiment performed toevaluate the use of a linear model.

FIG. 10 is a conceptual illustration of another example model used forthe quantitative analysis of biomedical image data.

FIGS. 11-12 show the results of additional example computer simulationsperformed to evaluate the use of a linear model.

FIG. 13 shows the results of an additional example experiment performedto evaluate the use of a linear model.

FIG. 14 is a flow chart of another example method that may be used withthe example method of FIG. 1.

FIG. 15 is a flow chart of still another example method that may be usedwith the example method of FIG. 1.

FIG. 16 is a simplified block diagram of an example system in which thepresent methods can be implemented.

FIG. 17 is a simplified block diagram of an example computing device.

FIG. 18 is a schematic illustrating a conceptual partial view of anexample computer program product that includes a computer program forexecuting a computer process on a computing device.

DETAILED DESCRIPTION

In the following detailed description, reference is made to theaccompanying figures, which form a part thereof. In the figures, similarsymbols typically identify similar components, unless context dictatesotherwise. The illustrative embodiments described in the detaileddescription, figures, and claims are not meant to be limiting. Otherembodiments may be utilized, and other changes may be made, withoutdeparting from the spirit or scope of the subject matter presentedherein. It will be readily understood that aspects of the presentdisclosure, as generally described herein, and illustrated in thefigures, can be arranged, substituted, combined, separated, and/ordesigned in a wide variety of different configurations, all of which areexplicitly contemplated herein.

Further, for illustration, certain aspects of the disclosure herein willbe described with respect to magnetic resonance imaging. It should beunderstood, however, that part or all of the described systems andmethods may apply equally to other types of biomedical imaging (e.g.,optical imaging, ultrasound imaging, etc.). Therefore, the describedembodiments should not be taken to be limiting.

I. EXAMPLE METHODS

Described herein are new methods and systems for analyzing models ofphenomena that change biomedical images. For purposes of example andexplanation, the methods and systems described herein may be used todetermine the value of at least one pharmacokinetic (PK) parameter thatis indicative of a perfusion, permeability, flow, extracellular volumefraction, or other physiological parameter of a tissue of interest(e.g., a tumor, muscle, organ, or other tissue). The determined value ofthe PK parameter may then be used to, for example, diagnose how thetissue of interest is responding to treatment, evaluate an abnormalcondition of blood perfusion, or prognosticate an outcome, among otheruses. Other aspects, uses, and advantageous of the methods and systemsdescribed herein are described further below.

With reference to FIG. 1, an example method 100 for determining a valueof a PK parameter is described. As shown in FIG. 1, initially at block102, the method 100 includes determine, based on a set of biomedicalimages, a first concentration-activity curve and a second-concentrationactivity curve. At block 104, the method 100 includes determine a valueof at least one pharmacokinetic (PK) parameter based on the firstconcentration-activity curve and the second concentration-activitycurve. Determining the value of the at least one PK parameter mayinclude determining a linear model that relates the firstconcentration-activity curve to the second concentration-activity curveat block 106, and at block 108, determining the value of the at leastone PK parameter based on application of a linear least squares (LLSQ)algorithm to the linear model. At block 110, the method 100 includescause a graphical display to provide a visual indication of the value ofthe at least one PK parameter. These steps are explained in thefollowing subsections.

Generally, the methods and functions described herein may be carried outby a computing system, such as computing system 1602 of FIG. 16described further below. Again, however, it should be understood thatthe computing system 1602 of FIG. 16 is set forth for purposes ofexample and explanation only, and should not be taken to be limiting.The present methods and functions may just as well be carried out inother systems having other arrangements.

a. Determine Concentration-Activity Curves

At block 102, the method 100 includes determine, based on a set ofbiomedical images, a first concentration-activity curve and a secondconcentration-activity curve. Each of the first concentration-activitycurve and the second concentration-activity curve may indicate aconcentration of a respective contrast agent within a respective regionof tissue.

In one example, the set of biomedical images may be a temporally dynamicseries of biomedical images. The set of biomedical images may includeany type of biomedical images or biomedical image data. For example, theset of biomedical images may be any type of MR images (e.g., F-MRI,DCE-MRI, ¹⁹F-DCE-MRI images), optical images, ultrasound images,positron emission tomography (PET) images, Computed Tomography images,or X-ray images, among other possibilities.

The respective contrast agent may be a substance used to enhance thecontrast of structures or fluids within a body for medical imaging. Theparticular type of contrast agent may vary depending on the type ofbiomedical imaging technique used to obtain the biomedical images. Inone case, the biomedical images may be MRI images, and the respectivecontrast agent may be a paramagnetic contrast agent. In another case,the biomedical images may be optical images, and the respective contrastagent may be a colored dye. In still another case, the biomedical imagesmay be ¹⁹F-MRI images, and the respective contrast agent may be ananoemulsion, such as a nanoemulsion that includes perfluorocarbons. Inyet another case, the contrast agent may be a substance that isresponsive to a biomarker. For instance, the contrast agent may beresponsive to zinc ions or enzyme activities. Other examples also exist.

A concentration-activity curve may generally be data that indicates theconcentration of a respective contrast agent as a function of timewithin a respective region of tissue (e.g., a voxel or group of voxels).As described further below, depending on the particular example, thefirst concentration-activity curve may indicate the concentration of thesame contrast agent or different respective contrast agents.Additionally, as described below, the concentration-activity curves mayindicate the concentration of a contrast agent within the same region oftissue or different respective regions of tissue.

Before turning to a more detailed description of the determination of aconcentration-activity curve, a brief overview of modeling thepharmacokinetics of a contrast agent is provided with reference to FIG.2. Such a model of the pharmacokinetics of a contrast agent mayultimately be used to generate a concentration-activity curve.

FIG. 2 is a conceptual illustration 200 of an example PK model used forthe quantitative analysis of biomedical image data. Fittingconcentration-curves to a PK model may be used to study the underlyingphysiology and/or pathology of an interest region of tissue, hereinreferred to as a tissue of interest (TOI). Parameters that describe themovement of a contrast agent across a vascular endothelium may include,for example: inflow rate from plasma to the TOI (K^(trans)), outflowrate from the TOI to plasma (k_(ep)), intravascular volume fraction(v_(p)), and extracellular-extravascular (EES) volume fraction (v_(e)).

PK models often make two assumptions. A first assumption is that a humanbody may be represented by one or more compartments, into and out ofwhich a contrast agent dynamically flows. A second assumption is thateach compartment is well-mixed such that a contrast agent entering acompartment is immediately distributed uniformly throughout thecompartment.

The model shown in FIG. 2 is a two-compartment model including a centralcompartment 202 and peripheral compartment 204. Intravascularextracellular volume fraction (plasma) is considered to be the centralcompartment 202 and EES volume fraction (tissue) is considered to theperipheral compartment 204.

According to the model, a contrast agent may be introduced into thevasculature and diffuses into the EES in a reversible processcharacterized by a distribution rate constant (K^(trans)) and aredistribution rate constant (k_(ep)=K^(trans)/v_(e)). Equation 1describes the process:

$\begin{matrix}{\frac{{{Ct}(t)}}{t} = {{K^{trans} \cdot {{Cp}(t)}} - {\frac{K^{trans}}{v_{e}} \cdot {{Ct}(t)}}}} & (1)\end{matrix}$

where Ct(t) is the concentration-activity curve for the tissue and Cp(t)is the concentration-activity curve in the plasma. The solution toEquation 1 is:

Ct(t)=K ^(trans)·∫₀ ^(T) Cp(t)·e ^((−K) ^(trans) ^(/v) ^(e) ^()·(T−t))dt  (2)

Equation 1 assumes that the tissue does not contain any vascular space,which may be correct for some tissues (e.g., muscle). However, in someexamples, ignoring the fraction of vascular space (v_(p)) may lead toerroneous results. Including v_(p) in the model increases accuracy andleads to the solution to Equation 1 shown by Equation 3:

Ct(t)=K ^(trans)·∫₀ ^(T) Cp(t)·e ^((−K) ^(trans) ^(/v) ^(e) ^()·(T−t))dt+v _(p)  (3)

In an example in which the biomedical images are DCE-MRI images,Equations 2 or 3 may be used to determine the firstconcentration-activity curve and the second concentration-activitycurve. The concentration-activity curve for plasma, Cp(t), is referredto as the arterial input function (AIF), and may be challenging todetermine Example approaches include taking arterial blood samplesduring acquisition of biomedical images or estimating the contrast agentconcentration from biomedical imaging data when a large vessel islocated in the same field of view of the analyzed tissue. A thirdexample approach is to eliminate the need of measuring the AIF by usinga second tissue, a reference region (RR), as a surrogate for the AIF.The third example approach, referred to as the RR model (RRM) is furtherdescribed below.

In an example in which the biomedical images are optical images, thefirst concentration-activity curve and the second concentration-activitycurve may be determined based on a pixel intensity value of a particularcolor or hue that corresponds to the color of a respective contrastagent. For instance, image processing may be used to determine theintensity of a particular color within a pixel or group of pixels withineach image, and the intensity may be converted to a concentration of thecontrast agent to determine a concentration-activity curve.

In an example in which a contrast agent is responsive to a biomarker,the first concentration-activity curve and the secondconcentration-activity curve may be determined based on changes inimaging signal, such as the MRI signal or the optical imaging signal,and whereby the biomarker-responsive agent changes the imaging signal ata different rate than the agent that is unresponsive to the biomarkerbut is responsive to the other conditions that also potentially affectthe first biomarker.

b. Determine a Value of at Least One PK Parameter

At block 104, the method 100 includes determine a value of at least onePK parameter based on the first concentration-activity curve and thesecond concentration-activity curve. According to the method 100,determining the value of the at least one PK parameter may includedetermining a linear model that relates the first concentration-activitycurve to the second concentration-activity curve at block 106, and atblock 108, determining the value of the at least one PK parameter basedon application of a linear least square fitting (LLSQ) algorithm to themodel.

As described previously, the PK parameter may be a parameter that isindicative of perfusion, permeability, flow, extracellular volumefraction, or other physiological parameter of a tissue of interest. Forinstance, the PK parameter may include a volume transfer constant, anextravascular extracellular volume fraction, or a redistribution rateconstant as described below, among other possibilities. Additionally,the PK parameter may be a parameter of a linear model that relates thefirst concentration-activity curve and the second concentration-activitycurve. The characteristics and derivations of example linear models,including a linear reference region model (LRRM) and a reference agentmodel (RAM), are now described in detail.

i. Linear Reference Region Model (LRRM)

The LRRM is one example of a linear model that may be determined atblock 106. LRRM is an improvement on the reference region model (RRM).Therefore, for purposes of example and explanation, a description of theRRM will first be described with reference to FIG. 3.

FIG. 3 is a conceptual illustration 300 of another example model thatmay be used in some examples for the quantitative analysis of biomedicalimage data. As mentioned briefly above, the example model shown in FIG.3 may be referred to as a reference region model (RRM), which eliminatesthe need of measuring the AIF by using a second tissue, a referenceregion (RR), as a surrogate for the AIF. According to the model, acontrast agent is introduced into the vasculature and diffuses into theEES of both a tissue of interest and a reference region of tissue in areversible process. The process can be described by Equations 4 and 5:

$\begin{matrix}{\frac{{C_{TOI}(t)}}{t} = {{K^{{trans},{TOI}} \cdot {{Cp}(t)}} - {\frac{K^{{trans},{TOI}}}{v_{e,{TOI}}} \cdot {C_{TOI}(t)}}}} & (4) \\{\frac{{C_{RR}(t)}}{t} = {{K^{{trans},{RR}} \cdot {{Cp}(t)}} - {\frac{K^{{trans},{RR}}}{v_{e,{RR}}} \cdot {C_{RR}(t)}}}} & (5)\end{matrix}$

Equations 4 and 5 can be combined to yield an equation that does notrequire Cp(t) as an input:

$\begin{matrix}{{C_{TOI}(t)} = {{\frac{K^{{trans},{TOI}}}{K^{{trans},{RR}}} \cdot {C_{RR}(t)}} + {\frac{K^{{trans},{TOI}}}{K^{{trans},{RR}}} \cdot \lbrack {k_{{ep},{RR}} - k_{{ep},{TOI}}} \rbrack \cdot {\int_{0}^{T}{{{C_{RR}(t)} \cdot ^{{({- k_{{ep},{TOI}}})} \cdot {({T - t})}}}\ {t}}}}}} & (6)\end{matrix}$

The physiological interpretation for K^(trans) varies depending on thebalance between capillary permeability (PS) and blood flow (F) in thetissue of interest. For high permeability situations (PS>>F):

K ^(trans) =F·ρ·(1−Hct)  (7)

where ρ is the tissue density and Hct is the hematocrit. When thetransport of the contrast agent is limited by permeability (PS<<F):

K ^(trans) =PS·ρ  (8)

In a third scenario, where contrast agent flux is limited by both flowand permeability:

K ^(trans)=(e ^(−PS/[F·(1−Hct)]) F·ρ·(1−Hct)  (9)

Using a RR model (RRM) slightly relaxes the requirements on temporalsampling for acquisition of biomedical images by using the concentrationof a contrast agent in a reference region as a surrogate for the AIF.Using the RRM, however, still involves previous knowledge or assumptionof the volume transfer constant and extravascular extracellular volume,K^(trans) and v_(e) respectively, of the reference region. In someexamples, such estimation may introduce systematic error into thedetermination of other PK parameters.

Another limitation of the RRM involves the use of a non-linear leastsquares (NLSQ) to estimate K^(trans) and v_(e) of the tissue of interestby fitting experimental concentration-activity curves to the RRM. Insome examples, NLSQ fittings perform poorly under low signal-to-noiseratio (SNR) conditions that are typical in clinical studies. Further,NLSQ fittings are highly dependent on the initial estimates of theconstants to be calculated. For example, for a heterogeneous tissue, theinitial estimates may not apply to all regions of the tissue, which maylead to inaccurate fittings. Additionally, the time required to obtain aparametric map through voxel-wise calculations that use a NLSQ algorithmcan often exceed one hour, which is impractical for most radiologyclinics.

In order to overcome the limitations described above, a linear RRM(LRRM) that may be used with the current method 100 was developed. TheLRRM is one example of a linear model that may be determined at block106.

The RRM may be formulated in terms of ratios of K^(trans,TOI) toK^(trans,RR) (known as R^(Ktrans)), K^(trans,TOI) to v_(e,TOI) (known ask_(ep,TOI)), and K^(trans,RR) to v_(e,RR) (known as k_(ep,RR)).Calculating R^(Ktrans), k_(ep,TOI), and k_(ep,RR) by applying a linearleast squares (LLSQ) fitting algorithm to a linear RRM removes the needfor an initial estimate of parameters. Using LLSQ can typically produceparametric maps within seconds, which can greatly facilitate thepractical use of these algorithms, for example, in clinical settings.Fitting concentration-activity curves to a LRRM also leads to betterperformance compared to a using a NLSQ algorithm to fitconcentration-activity curves to a non-linear RRM (NRRM). Additionally,using the LRRM described herein provides for more precise and accurateresults under conditions of low SNR and/or slow temporal resolution thana NRRM.

Given the behavior of a contrast agent within a tissue of interest and areference region, described by Equations 10 and 11:

$\begin{matrix}{\frac{{C_{TOI}(t)}}{t} = {{K^{{trans},{TOI}} \cdot {{Cp}(t)}} - {K_{{ep},{TOI}} \cdot {C_{TOI}(t)}}}} & (10) \\{\frac{{C_{RR}(t)}}{t} = {{K^{{trans},{RR}} \cdot {{Cp}(t)}} - {K_{{ep},{RR}} \cdot {C_{RR}(t)}}}} & (11)\end{matrix}$

it is possible to solve equation 11 for Cp(t), yielding

$\begin{matrix}{{{Cp}(t)} = {{\frac{1}{K^{{trans},{RR}}} \cdot \frac{{C_{RR}(t)}}{t}} + {\frac{1}{v_{e,{RR}}} \cdot {C_{RR}(t)}}}} & (12)\end{matrix}$

Substituting Equation 12 into Equation 10 eliminates the dependence ofEquation 10 on Cp(t) and yields Equation 13:

$\begin{matrix}{\frac{{C_{TOI}(t)}}{t} = {K^{{trans},{TOI}} \cdot \lbrack {{\frac{1}{K^{{trans},{RR}}} \cdot \frac{{C_{RR}(t)}}{t}} + {\frac{1}{v_{e,{RR}}} \cdot {C_{RR}(t\rbrack}} - {k_{{ep},{TOI}} \cdot {C_{TOI}(t)}}} }} & (13)\end{matrix}$

As previously mentioned, the solution to Equation 13 is:

$\begin{matrix}{{C_{TOI}(t)} = {{\frac{K^{{trans},{TOI}}}{K^{{trans},{RR}}} \cdot {C_{RR}(t)}} + {\frac{K^{{trans},{TOI}}}{K^{{trans},{RR}}} \cdot \lbrack {k_{{ep},{RR}} - k_{{ep},{TOI}}} \rbrack \cdot {\int_{0}^{T}{{{C_{RR}(t)} \cdot ^{{({- k_{{ep},{TOI}}})} \cdot {({T - t})}}}\ {t}}}}}} & (14)\end{matrix}$

and can be used to estimate the relative transfer constant of TOIcompared to RR, R^(Ktrans).

A LRRM can be derived by integrating both sides of Equation 13 andassuming that the initial concentration of the contrast agent in bothtissues is equal to zero to obtain Equation 15:

$\begin{matrix}{{C_{TOI}(t)} = {{R^{Ktrans} \cdot {C_{RR}(t)}} + {\frac{K^{{trans},{TOI}}}{v_{e,{RR}}} \cdot {\int_{0}^{T}{{C_{RR}(t)}\ {t}}}} - {k_{{ep},{TOI}} \cdot {\int_{0}^{T}{{C_{TOI}(t)}\ {t}}}}}} & (15)\end{matrix}$

Expressed in matrix form, Equation 14 is equal to:

$\begin{matrix}{\overset{harpoonup}{A} = {\overset{harpoonup}{M} \cdot \overset{harpoonup}{b}}} & (16) \\{where} & \; \\{{{{{\overset{harpoonup}{A} = \begin{pmatrix}{C_{TOI}( {T\; 1} )} \\{C_{TOI}( {T\; 2} )} \\\vdots \\{C_{TOI}( {T\; n} )}\end{pmatrix}};}\overset{harpoonup}{M} = \begin{pmatrix}{C_{RR}( {T\; 1} )} & {\int_{0}^{T\; 1}{{C_{RR}(t)}\ {t}}} & {- {\int_{0}^{T\; 1}{{C_{TOI}(t)}\ {t}}}} \\{C_{RR}( {T\; 2} )} & {\int_{0}^{T\; 2}{{C_{RR}(t)}\ {t}}} & {- {\int_{0}^{T\; 2}{{C_{TOI}(t)}\ {t}}}} \\\vdots & \vdots & \vdots \\{C_{RR}( {T\; n} )} & {\int_{0}^{Tn}{{C_{RR}(t)}\ {t}}} & {- {\int_{0}^{T\; n}{{C_{TOI}(t)}\ {t}}}}\end{pmatrix}};}{\overset{harpoonup}{b} = ( \frac{\begin{matrix}R^{Ktrans} \\K^{{trans},{TOI}}\end{matrix}}{\begin{matrix}v_{e,{RR}} \\k_{{ep},{TOI}}\end{matrix}} )}} & (17)\end{matrix}$

If the elements of columns two and three of vector M are approximated bynumerical integration, Equation 16 can be solved for the elements of busing conventional LLSQ. In one example, the application of a LLSQalgorithm of block 108, may include applying a LLSQ algorithm to theLRRM model described by Equations 15 and 16. For instance, applying theLLSQ to Equations 15 and 16 may yield the PK parameters of vector b ofEquations 15 and 16.

R^(Ktrans), K^(trans,TOI), K^(trans,RR), k_(ep,TOI), and v_(e,TOI),v_(e,RR) are a examples of PK parameters that may be determined based ona linear model, the LRRM. In an example in which the LRRM is used as thelinear model of the method 100, C_(TOI)(t) and C_(RR)(t) may be thefirst concentration-activity curve and the second concentration-activitycurve. Thus, the first concentration-activity curve may indicate theconcentration of a first contrast agent within an interest region oftissue, the TOI, and the second concentration-activity curve mayindicate the concentration of the first contrast agent within areference region of tissue, the RR.

ii. Simulation and Experimental Results Using LRRM

Computer simulations and experiments were performed to test the use ofthe LRRM in accordance with the method 100. Advantageously, the resultsindicated that using the LRRM is faster than the NRRM and relaxesrequirements for temporal resolution and SNR. The simulation results andexperimental results for the LRRM are described below. The simulationand experimental results are provided as examples. Other results arealso possible and the results presented herein should not be taken to belimiting in any way.

Computer simulations were used to study the effect of injection speed,temporal resolution, and statistical noise on the accuracy and precisionof LRRM and NRRM to estimate R^(Ktrans) and k_(ep,TOI). Three arterialinput functions (AIFs) with three different injection speeds weresimulated using Equation 18

Cp(t)=A·t ^(B) ·e ^((−t·C)) +D·[1−e ^((−t·E)) ]·e ^((−t·F))  (18)

and the following parameters for each speed: A) Fast 10-sec. injectionspeed, A=30.0 mM, B=1.0, C=4.00 min⁻¹, D=0.65 mM, E=5.0 min⁻¹; F=0.04min⁻¹; B) Moderate 40-sec. injection speed, A=120.0 mM, B=3.0, C=4.34min⁻¹, D=0.80 mM, E=1.0 min⁻¹, F=0.07 min⁻¹; C) Slow 80-sec. injectionspeed, A=55.0 mM, B=5.5, C=3.98 min⁻¹, D=0.80 mM, E=1.0 min⁻¹, F=0.07min⁻¹ Each of the AIFs were calculated from t=0 to t=35 minutes at atemporal resolution of 1.0 seconds. A series of C_(TOI)(T) and C_(RR)(T)curves were then calculated using each AIF and Equation 19:

C _(t)(T)=K ^(trans)∫₀ ^(T) C _(p)(t)·e ^([−k) ^(ep) ^(·(T−t)]) dt  (19)

with K^(trans)=0.25 min⁻¹ and k_(ep)=0.62 min⁻¹ for the TOI curve andK^(trans)=0.10 min⁻¹ and k_(ep)=1.00 min⁻¹ for the RR curve.

To study systematic error, all the AIFs, C_(TOI)(T), and C_(RR)(T) weredown sampled from their initial resolution of 1 second to a newresolution of 60 seconds in intervals of 1 second. Further, the AIF withthe best performance in the systematic error analysis was used to studythe combined effect of SNR and temporal resolution on the accuracy andprecision of the LRRM and NRRM. To do so, all curves were recalculatedat Δt=1, 10, 30, and 60 seconds. Also, Gaussian noise was added toproduce SNRs from 5 to 50. Each simulation was then repeated 1000 timesfor each condition of SNR and Δt. The accuracy of parameters wasevaluated in terms of the mean percentage error and the precision wasevaluated in terms of the standard deviation of the percentage error.

The estimated values of R^(Ktrans) and k_(ep,TOI) from these simulationswere compared to their initial values used to create the simulated data.FIG. 4 shows the systematic errors as a function of temporal resolutionand injection speed. As shown in FIG. 4, the NRRM systematicallyunderestimated R^(Ktrans) and overestimated k_(ep,TOI), and thesesystematic errors became worse as the temporal sampling became slower.The results shown in FIG. 4 indicate that DCE-MRI with LRRM can produceaccurate results with a moderate or fast injection speed, while DCE-MRIwith NRRM cannot produce accurate results unless the temporal samplingis extremely fast.

FIG. 5 shows the error in estimated R^(Ktrans) and k_(ep,TOI) as afunction of SNR. As shown in FIG. 5, both the LRRM and NRRM showedgreater average errors in estimated R^(Ktrans) and k_(ep,TOI) as the SNRdecreased. Yet, the errors with LRRM were smaller than the errors withNRRM for all SNRs. The standard deviations of these results were alsosmaller when using LRRM relative to NRRM at SNR>10. Also, for a SNR of60 these standard deviations were negligible for NRRM while thesestandard deviations were substantial for NRRM, indicating the LRRMprovided more precise fittings to the data.

FIG. 6 shows the error in estimated R^(Ktrans) as a function of SNR andtemporal resolution. As shown in FIG. 6, the errors in estimatingR^(Ktrans) increased as the SNR decreased for both LRRM and NRRM. Yet,the errors in estimated R^(Ktrans) with LRRM were smaller than theerrors with NRRM for all SNRs. Also, the errors in estimated R^(Ktrans)with NRRM were strongly dependent on temporal resolution while estimateswith LRRM were relatively insensitive to temporal resolution. Theresults shown in FIG. 6 illustrate that DCE-MRI with LRRM can producemore precise results for estimated R^(Ktrans) than NRRM, especially atlow SNRs, and regardless of injection speeds.

FIG. 7 shows the error in estimated k_(ep,TOI) as a function of SNR andtemporal resolution. As shown in FIG. 7, the errors in estimatedk_(ep,TOI) with LRRM were also smaller than errors with NRRM for fasttemporal resolutions of 5 and 10 seconds, but the errors with LRRM werelarger than errors with NRRM for slower temporal resolutions. Both theLRRM and NRRM produced errors in estimated k_(ep,TOI) that weredependent on temporal resolution. These results demonstrated thatDCE-MRI with LRRM can produce more precise results for estimatedk_(ep,TOI) than the NRRM at fast temporal resolutions.

The results of experiments also indicated that the LRRM may be moresuitable than the NRRM for many clinical studies. For instance, resultsindicated that NRRM may only produce accurate and precise results whenthe temporal resolution is fast, which may not be feasible in allradiology clinics.

To experimentally compare the results of LRRM and NRRM, a female severecombined immunodeficiency (SCID) mouse weighing 21 grams was housed andmaintained under specific pathogen-free conditions. The mouse wasinjected subcutaneously in the right flank with 10×10⁶ MDA-231 cells in0.1 mL saline. The tumor was measured every 2 to 3 days using electroniccalipers and allowed to grow for 22 days to an average volume of 250 mm³before initiating MRI studies.

The mouse was scanned 22 days after tumor implantation using a BrukerBiospec 7.0 T scanner that was equipped with a 72 mm birdcage coil(Bruker Biospin, Inc.). Before the imaging scan, the mouse wasanesthetized with 1.5-2% isoflurane in O₂ carrier gas, a 27 G catheterwas inserted in the tail vein and physiological monitoring leads wereconnected to monitor respiration rate and core body temperature duringthe MRI session. The mouse was kept at 37.0±0.2° C. during the MRIstudies using warmed air that was controlled by an automatedtemperature-feedback system (SA Instruments, Inc.).

DCE-MRI studies were conducted by acquiring a set of 3 contiguous axialimage slices through the xenograft flank tumor, and 3 axial image slicesthrough the femoral artery. The orientation of these image slices weredetermined by acquiring T₂-weighted (T_(2w)) MR images with a RARE MRIprotocol that used the following parameters: TR=1.0 sec, TE=8.2 msec,slice thickness=0.5 mm; matrix=128×128, FOV=3.5 cm², in-planeresolution=273 μm²; NEX=1, and RARE=4.

A parametric map of endogenous T₁ relaxation time, R₁(0) with the samegeometry as the T_(2w) images, was obtained by acquiring a series ofRARE-variable TR (RARE-VTR) images with the following parameters:TR=0.375, 0.75, 1.5, 3.0, and 6.0 sec, TE=9.07 msec, NEX=1, RAREfactor=2, and fitting the data to Equation 21:

$\begin{matrix}{S = {A \cdot \lbrack \frac{\sin \; {\alpha ( {1 - ^{({{- {TR}} \cdot R_{1}})}} )}}{1 - {{^{({{- {TR}} \cdot R_{1}})} \cdot \cos}\; \alpha}} \rbrack}} & (21)\end{matrix}$

where α is the flip angle (90 degrees in the experiment), and A is aconstant that includes the proton density and scanner gain. It wasassumed that TE<<T₂*.

A series of T₁-weighted (T_(1w)) images were acquired for DCE-MRIanalysis using a RARE MRI protocol with the same geometry as the R₁(0)maps, and the following parameters was used: TR=250 msec, TE=8.2 msec,NEX=2 averages, a RARE factor=2. A total of 65 image sets were acquiredfor a total acquisition time of 34.66 min and a temporal resolution of32 sec/image. After the fifth image set was acquired, 0.2 mmol/KgGd-DTPA (Magnevist™, Bayer Healthcare) was injected through the tailvein catheter in a total volume of 0.25 mL during 60 seconds. Finally,the R₁(0) maps and Equation 21 were used to estimate ΔR₁(t) for eachvoxel within the TOI and the average ΔR₁(t) of the RR.

$\begin{matrix}{{R_{1}(t)} = {\frac{1}{- {TR}} \cdot {\ln \lbrack {1 - {\frac{S(t)}{S(0)}( {1 - ^{({{- {TR}} \cdot {R_{1}{(0)}}})}} )}} \rbrack}}} & (22)\end{matrix}$

The ΔR₁(t) curves for the RR and TOI were fitted to the LRRM and NRRM togenerate tumor parametric maps of R^(Ktrans) and k_(ep,TOI) Theconcentration of a contrast agent may be determined based on the effectof the contrast agent on water relaxation. The T₁ relaxation time andthe concentration of the contrast agent are related through thefollowing equations:

$\begin{matrix}{{R_{1}(t)} = {{r \cdot {{CA}(t)}} + {R_{1}(0)}}} & (23) \\{\frac{{R_{1}(t)} - {R_{1}(0)}}{r} = {\frac{\Delta \; {R_{1}(t)}}{r} = {{CA}(t)}}} & (24)\end{matrix}$

where R₁(t) is the transverse relaxation rate (1/T₁) at each time pointof the DCE-MRI experiment, CA(t) is the contrast agent concentration atthat same time point, R₁(0) is the native R₁ of the tissue beforeinjection of the contrast agent, and r is the contrast agent relaxivity.

The time to calculate a parametric map of R^(Ktrans) for a 128×128 imagematrix was approximately 6.64 seconds using LRRM and 68.88 minutes usingNRRM, using a typical notebook computer. Calculating a parametric map ofk_(ep,TOI) required similar calculation times. Therefore, the LRRMproduced results approximately 620 times faster than NRRM.

FIG. 8 shows an experimental evaluation of the LRRM and NRRM. Using boththe LRRM and NRRM to solve for R^(Ktrans) and k_(ep,TOI) and k_(ep,RR)yielded results having a root mean squared error (RMSE) of 0.008 mM. Theestimated parameters were comparable to those reported for other tumortypes and muscle. As shown in FIG. 8, the LRRM and NRRM performedsimilarly when fitting data with high SNR, moderate temporal resolution,and moderate injection speed.

FIG. 9 shows the effect of temporal resolution on the accuracy ofestimated R^(Ktrans) and k_(ep,TOI). The ΔR₁ curves for each voxel andthe reference region were down-sampled from their initial values of 32sec to values of 64, 96, and 128 sec. R^(Ktrans) and k_(ep,TOI) wereestimated with LRRM and NRRM at each temporal resolution. In FIG. 9,scatter plots showing the values at estimated at the slower temporalresolutions compared with values at 32 sec are presented. As shown inFIG. 9, when LRRM was used, R^(Ktrans) and k_(ep,TOI) estimated withslower resolution were highly correlated with R^(Ktran,32s) andk_(ep,TOI,32), while the NRRM estimated different R^(Ktrans) andk_(ep,TOI) at each temporal resolution. The trends are furtherillustrated with reference to the Spearmen correlation coefficients forthe scatter plots that are shown in FIG. 9. For example, R^(Ktrans,32)was highly correlated with R^(Ktrans,128) when the LRRM was used(r=0.91), while the R^(Ktrans,32) and R^(Ktrans,128) estimated with NRRMshowed a lower correlation (r=0.61). Thus, the experimental results showthat LRRM is relatively insensitive to temporal resolution, while NRRMis more sensitive to temporal resolution.

iii. Reference Agent Model (RAM)

Another example of a linear model that may be determined at block 106 isthe reference agent model. A detailed description of the RAM is nowdescribed.

In some applications, such as DCE-MRI, for example, the inability toaccount for an exact value of hematocrit in blood may lead to improperlyestimated values of Cp(t) and K^(trans). For example, most DCE-MRIstudies assume that the arterial hematocrit is 40%. However, thehematocrit decreases as vessel diameter decreases, so that hematocrit intumor microvasculature, for example, ranges between 20%-80%.

The RAM compares the pharmacokinetics of two contrast agents in the sametissue, so that one contrast agent may be used as a reference for thesecond contrast agent. Because both contrast agents are located in thesame tissue and experience the same hematocrit, a ratiometric comparisonof biomedical images of both contrast agents is independent of thehematocrit. More generally, a ratiometric approach has potential tocancel other physiological characteristics as well.

The differential equations that describe the pharmacokinetic behavior oftwo contrast agents within the same tissue are described with respect toFIG. 10. FIG. 10 is a conceptual illustration 1000 of another examplemodel used for the quantitative analysis of biomedical image data.Similar to the conceptual illustration 300 and example model of FIG. 3,two assumptions are made. A first assumption is that a human body may berepresented by one or more compartments, into and out of which acontrast agent dynamically flows. A second assumption is that eachcompartment is well-mixed such that a contrast agent entering acompartment is immediately distributed uniformly throughout thecompartment.

The model shown in FIG. 10 is a two-compartment model including acentral compartment 1002 and peripheral compartment 1004. Intravascularextracellular volume fraction (plasma) is considered to be the centralcompartment 1002 and EES volume fraction (tissue) is considered to bethe peripheral compartment 1004.

According to the RAM, a first contrast agent (CA-1) and a secondcontrast agent (CA-2) are introduced into the vasculature and diffuseinto the EES in a reversible process characterized by a distributionrate constant (K^(trans)) and a redistribution rate constant(k_(ep)=K^(trans)/v_(e)). The following differential equations describethe process:

$\begin{matrix}{\frac{{{Ct}_{{CA} - 1}(t)}}{t} = {{K^{{trans},{{CA} - 1}} \cdot \frac{{Cb}_{{CA} - 1}(t)}{( {1 - {Hct}} )}} - {\frac{K^{{trans},{{CA} - 1}}}{v_{e}} \cdot {{Ct}_{{CA} - 1}(t)}}}} & (25) \\{\frac{{{Ct}_{{CA} - 2}(t)}}{t} = {{K^{{trans},{{CA} - 2}} \cdot \frac{{Cb}_{{CA} - 2}(t)}{( {1 - {Hct}} )}} - {\frac{K^{{trans},{{CA} - 2}}}{v_{e}} \cdot {{Ct}_{{CA} - 2}(t)}}}} & (26)\end{matrix}$

Equation 27 relates the concentration of a contrast agent in blood tothe concentration of the contrast agent in plasma.

$\begin{matrix}{{{Cp}(t)} = \frac{{Cb}(t)}{( {1 - {Hct}} )}} & (27)\end{matrix}$

Assuming that Cb_(CA-1)(t) is equal to Cb_(CA-2)(t), Equation 25 can besolved for Cp(t):

$\begin{matrix}{\frac{{Cb}(t)}{( {1 - {Hct}} )} = {{\frac{1}{K^{{trans},{{CA} - 2}}} \cdot \frac{{{Ct}_{{CA} - 2}(t)}}{t}} + {\frac{1}{v_{e}} \cdot {{Ct}_{{CA} - 2}(t)}}}} & (28)\end{matrix}$

Substituting Equation 27 into Equation 24 then eliminates the dependenceof the analysis on Cb(t) and Hct:

$\begin{matrix}{\frac{{{Ct}_{{CA} - 1}(t)}}{t} = {{\frac{K^{{trans},{{CA} - 1}}}{K^{{trans},{{CA} - 2}}} \cdot \frac{{{Ct}_{{CA} - 2}(t)}}{t}} + {\frac{K^{{trans},{{CA} - 1}}}{v_{e}} \cdot \lbrack {{{Ct}_{{CA} - 2}(t)} - {{Ct}_{{CA} - 1}(t)}} \rbrack}}} & (29)\end{matrix}$

Finally, integrating both sides of Equation 29 and assuming that theinitial concentrations of CA-1 and CA-2 are equal to zero yields theworking equation for the RAM:

$\begin{matrix}{{{Ct}_{{CA} - 1}(T)} = {{\frac{K^{{trans},{{CA} - 1}}}{K^{{trans},{{CA} - 2}}} \cdot {{Ct}_{{CA} - 2}(T)}} + {\frac{K^{{trans},{{CA} - 1}}}{v_{e}} \cdot \lbrack {{\int_{0}^{T}{{{Ct}_{{CA} - 2}(t)}\ {t}}} - {\int_{0}^{T}{{{Ct}_{{CA} - 1}(t)}\ {t}}}} \rbrack}}} & (30) \\{{{Ct}_{{CA} - 1}(T)} = {{R^{ktrans} \cdot {{Ct}_{{CA} - 2}(T)}} + {k_{{ep},{{CA} - 1}} \cdot \lbrack {{\int_{0}^{T}{{{Ct}_{{CA} - 2}(t)}\ {t}}} - {\int_{0}^{T}{{{Ct}_{{CA} - 1}(t)}\ {t}}}} \rbrack}}} & (31)\end{matrix}$

Expressed in matrix form, Equation 31 is equal to Equation 32:

$\begin{matrix}{\mspace{79mu} {\overset{harpoonup}{A} = {\overset{harpoonup}{M} \cdot \overset{harpoonup}{b}}}} & (32) \\{\mspace{79mu} {where}} & \; \\{\mspace{79mu} {{{\overset{harpoonup}{A} = \begin{pmatrix}{{Ct}_{{CA} - 1}( {T\; 1} )} \\{{Ct}_{{CA} - 1}( {T\; 2} )} \\\vdots \\{{Ct}_{{CA} - 1}( {T\; n} )}\end{pmatrix}};{\overset{harpoonup}{b} = \begin{pmatrix}R^{Ktrans} \\k_{{ep},{{CA} - 1}}\end{pmatrix}}}{{\overset{harpoonup}{M} = \begin{pmatrix}{{Ct}_{{CA} - 2}( {T\; 1} )} & \lbrack {{\int_{0}^{T\; 1}{{{Ct}_{{CA} - 2}(t)}\ {t}}} - {\int_{0}^{T\; 1}{{{Ct}_{{CA} - 1}(t)}\ {t}}}} \rbrack \\{{Ct}_{{CA} - 2}( {T\; 2} )} & \lbrack {{\int_{0}^{T\; 2}{{{Ct}_{{CA} - 2}(t)}\ {t}}} - {\int_{0}^{T\; 2}{{{Ct}_{{CA} - 1}(t)}\ {t}}}} \rbrack \\\vdots & \vdots \\{{Ct}_{{CA} - 2}( {T\; n} )} & \lbrack {{\int_{0}^{Tn}{{{Ct}_{{CA} - 2}(t)}\ {t}}} - {\int_{0}^{T\; n}{{{Ct}_{{CA} - 1}(t)}\ {t}}}} \rbrack\end{pmatrix}};}}} & (33)\end{matrix}$

The elements of column two of the vector M in Equation 32 can beapproximated by numerical integration. Equations 31 and 32, therefore,represent a system of linear equations, which can be solved for theelements of vector b. The RAM can also be expressed as a linearequation:

$\begin{matrix}{y = {{m \cdot x} + b}} & (34) \\{where} & \; \\{{{y = \frac{{Ct}_{{CA} - 1}(T)}{{Ct}_{{CA} - 2}(T)}};{x = \frac{\lbrack {{\int_{0}^{T}{{{Ct}_{{CA} - 2}(t)}\ {t}}} - {\int_{0}^{T}{{{Ct}_{{CA} - 1}(t)}\ {t}}}} \rbrack}{{Ct}_{{CA} - 2}(T)}};}{{m = k_{{ep},{{CA} - 1}}};{b = R^{Ktrans}}}} & (35)\end{matrix}$

Equations 34 and 35 can then be used to estimate R^(Krans) andk_(ep,CA-1) by linear regression fitting.

In one example, the application of a LLSQ algorithm to the linear modelin block 108 of the method 100 may include applying a LLSQ algorithm toEquations 31 and 32. In another example, the application of a LLSQalgorithm to the linear model in block 108 of the method 100 may includeapplying a LLSQ algorithm to Equations 34 and 35.

As mentioned briefly above, if the flux of a contrast agent across atissue endothelium has low permeability relative to flow, leans can beapproximated by Equation 8:

K ^(trans) =PS·ρ  (8)

where P is the total permeability of the capillary wall, S is thesurface area per unit mass of tissue, and ρ is the density of tissue.

In an example in which a biomedical imaging study is conducted withlarge molecule contrast agents (e.g., nanoemulsions of perfluorocarbons)that match the permeability-limited model, the ratio of K^(trans) valuesfor two large molecule contrast agents simultaneously detected in thesame tissue can be expressed as R^(Ktrans):

$\begin{matrix}{R^{Ktrans} = {\frac{K^{{trans},{{CA} - 1}}}{K^{{trans},{{CA} - 2}}} = {\frac{P_{{CA} - 1} \cdot S \cdot \rho}{P_{{CA} - 2} \cdot S \cdot \rho} = \frac{P_{{CA} - 1}}{P_{{CA} - 2}}}}} & (36)\end{matrix}$

where P_(CA-1) and P_(CA-2) are the apparent permeabilities of thetissue for CA-1 and CA-2 respectively.

R^(Ktrans), K^(trans,CA-1), K^(trans,CA-2), k_(ep,CA-1), k_(ep,CA-2),and v_(e) are all examples of PK parameters that may be determined basedon a linear model, the RAM. In an example in which the RAM is used asthe linear model of the method 100, Ct_(CA-1)(t) and Ct_(CA-2)(t) may bethe first concentration-activity curve and the secondconcentration-activity curve. Thus, the first concentration-activitycurve may indicate the concentration of a first contrast agent within aninterest region of tissue, the TOI, and the secondconcentration-activity curve may indicate the concentration of a secondcontrast agent within the interest region of tissue.

iv. Simulation and Experimental Results Using RAM

Computer simulations and experiments were performed to test the use ofthe RAM in accordance with the method 100. The simulation results andexperimental results for the RAM are described below. Again, thesimulation and experimental results are provided as examples. Otherresults are also possible and the results presented herein should not betaken to be limiting in any way.

Computer simulations were used to study the effect of temporalresolution and statistical noise on the accuracy and precision of RAM.The AIF for ¹⁹F nanoemulsions used in the study were experimentallymeasured and fitted to the following equation:

C _(p)(t)=A·e ^((−α·T)) +B·e ^((−β·T))  (36)

The estimated parameters α and β were combined with previously reportedvalues for A and B for nanoparticles of similar size to simulate the AIFfrom 0 to 30 min. The values for each parameter were A=24.0 KgL⁻¹;α=16.54 min⁻¹; B=5.5 KgL⁻¹; and β=0.03935 min⁻¹ Ct_(CA-1)(t) andCt_(CA-2)(t) curves were calculated using Equation 37:

$\begin{matrix}{{C_{t}(T)} = {K^{trans}{\int_{0}^{T}{{{C_{p}(t)} \cdot ^{\lbrack{\frac{K^{trans}}{v_{e}}{({T - t})}}\rbrack}}\ {t}}}}} & (37)\end{matrix}$

with K^(trans) equal to 8.2×10⁻³ min⁻¹ and 4.3×10⁻³ min⁻¹ for CA-1 andCA-2 respectively. The parameter v_(e) was set to 0.5 mL/mL for bothcurves.

To test the effect of temporal resolution, the AIF, Ct_(CA-1)(t), andCt_(CA-2)(t) were down-sampled from Δt=one second to Δt=0.5, 2.0, and4.0 min. The combined effect of SNR and Δt on accuracy and precision wasstudied by adding Gaussian noise to produce SNRs from 5 to 100.

FIG. 11 shows the systematic error in R^(Ktrans) and k_(ep,CA-1) as afunction of R^(Ktrans) and temporal resolution. As shown in FIG. 11, thesystematic error is greatest for slow temporal resolutions whenR^(Ktrans) is large. Additionally, when the PK of both contrast agentsis comparable, creating a small to moderate value of R^(Ktrans), a slowtemporal resolution can accurately R^(Ktrans) with a systematic error ofless than 1%. FIG. 11 also shows that the systematic error in estimatingk_(ep,CA-1) is greater than 1% for all cases except the fastest temporalresolution.

FIG. 12 shows the systematic error in R^(Ktrans) and k_(ep,CA-1) as afunction of SNR. As shown in FIG. 12, error in estimating R^(Ktrans) isless than 2% with SNR greater than 30:1, and less than 10% with SNRgreater than 15:1. Additionally, the accuracy and precision ofk_(ep,CA-1) is also dependent on SNR. However, errors in k_(ep,CA-1) aregreater than errors in R^(Ktrans), suggesting that estimating theratiometric parameter R^(Ktrans) has advantages relative to estimatingk_(ep,CA-1).

The results of experiments also indicated that RAM may be used for manyclinical studies. To experimentally evaluate the use of RAM, a mousemodel of a subcutaneous MDA-MB-231 tumor was used to perform in vivo¹⁹F-DCE-MRI. A total of 18 multi-echo images were acquired by iteratingthe selective detection of perfluorinated-15-crown-ether (PCE) andperfluorooctane (PFO), which produced 9 images of each nanoemulsion thatwere unshuffled to show the temporal dependence of the ¹⁹F MRI signalsfrom each nanoemulsion. The mixture of the nanoemulsions was injectedduring the acquisition of the second image, so that the first two imagesshowed no significant ¹⁹F MR signal. Subsequent images showed aconsistently increasing ¹⁹F MRI signal in the tumor. ¹⁹F MRI signals inthe tumor were converted to ¹⁹F concentrations, and the RAM was used toestimate a value of R^(Ktrans).

FIG. 13 illustrates example concentration-activity curves for PCE andPFO. The example concentration-activity curves indicate the ¹⁹Fconcentrations of the two nanoemulsions. As shown in FIG. 13, anR^(Ktrans) value of approximately 1 was determined. The R^(Ktrans) valuewas expected because of the similar relative diameters of PCE and PFO.The RMSE of the RAM fitting was 2.7×10⁻² mM, which is approximately 1%error relative to the magnitudes of the ¹⁹F concentrations that weremeasured. Additionally, k_(ep,CA-1) was estimated to be 9.1×10⁻³ min⁻¹,which was in the range of expected values. Therefore, the in vivo¹⁹F-DCE-MRI results demonstrated that two or more MRI contrast agentscan selectively be detected during a single MRI scan session andR^(Ktrans) can be measured with RAM.

c. Cause a Graphical Display to Provide a Value

At block 110, the method 100 includes cause a graphical display toprovide a visual indication of the value of the at least one PKparameter. The graphical display may be any type of display device, suchas a display device that is communicatively coupled to a computingsystem that is configured to perform the method 100. In one example, themethod 100 may include causing the graphical display to provide thevisual indication in response to the determining of the value of the PKparameter.

II. ADDITIONAL FUNCTIONS AND EXAMPLE APPLICATIONS

In some examples, the method 100 of FIG. 1 may further includeadditional blocks. FIGS. 14 and 15 are flow charts of example methods200 and 300. One or more blocks of the method 100 may be combined withany of the blocks of methods 1400 and/or 1500. For instance, the method1400 or the method 1500 may be performed after block 110 of FIG. 1.Other configurations are also possible. Alternatively, either of themethods 1400 and 1500 may be performed independently of the method 100.

As shown in FIG. 2, initially, at block 1402, the method 1400 mayinclude compare a determined value of at least one PK parameter to athreshold or a previous value. The determined value may be a valuedetermined using the method 100, for example, or may be a valuedetermined using other methods. In one example, the determined value maybe a volume transfer constant, and the volume transfer constant may becompared to a volume transfer constant threshold or to a prior value ofa volume transfer constant determined at an early time period. Thecomparison may determine, for instance, whether the determined value isgreater than or less than the threshold or previous value.

At block 1404, the method 1400 may include determine aresponse-indication corresponding to a treatment of an interest regionof tissue based on the comparison. For instance, the interest region oftissue may be a tumor, and the response-indication may be an indicationof whether the tumor is responding to an anti-angiogenic therapy. In oneexample, if the comparison indicates that a determined volume transferconstant is less than a previously determined volume transfer constant,the response-indication may be that the tumor is responding positivelyto the anti-angiogenic therapy.

The methods of the invention may thus be used, for example, to diagnosecancer in a subject. In one embodiment, the cancer comprises a solidtumor. Exemplary tumor types that can be diagnosed using the methods ofthe invention include, but are not limited to, lung cancer, breastcancer, bladder cancer, thyroid cancer, liver cancer, pleural cancer,pancreatic cancer, ovarian cancer, cervical cancer, testicular cancer,colon cancer, anal cancer, bile duct cancer, gastrointestinal carcinoidtumors, esophageal cancer, gall bladder cancer, rectal cancer, appendixcancer, small intestine cancer, stomach (gastric) cancer, renal cancer,cancer of the central nervous system, skin cancer, choriocarcinomas;head and neck cancers, osteogenic sarcomas, B-cell lymphoma,non-Hodgkin's lymphoma, Burkitt's lymphoma, fibrosarcoma, neuroblastoma,glioma, and melanoma. The subject may be any suitable subject, includinga mammal such as a human. The subject may be any subject with one ormore risk factors for cancer, or who is otherwise identified by anattending physician as in need of a test for diagnosing cancer. In oneexemplary embodiment, a subject “at risk of breast cancer” is anysubject considered to be in a risk group for breast cancer. In oneembodiment, the subject is a woman. In other embodiments, the subjecthas one or more of a lump in their breast tissue, lymph nodes, orarmpit; changes in breast size or shape; skin dimpling; nippleinversion; spontaneous single-nipple discharge; a family/personalhistory of breast cancer; or is a carrier of a mutation in the BRCA orother gene that predisposes one to breast cancer. Based on the teachingsherein and the level of skill in the art, it will be apparent to thoseof skill in the art how to identify subjects at risk of other cancers.

The methods of the invention may further be used, for example, to assessa response of a subject with cancer to a prescribed anti-cancertreatment, including but not limited to radiation therapy, chemotherapy(i.e., mechlorethamine, cyclophosphamide, chlorambucil, ifosfamide,azathioprine, mercaptopurine, taxanes, vincristine, vinblastine,paclitaxel, imatinib, gefitinib, etc.), anti-angiogenic therapy (ex:bevacizumab), hormonal therapy, monoclonal antibody therapy(trastuzumab, rituximab, etc.), and photodynamic therapy. In theseembodiments, the methods are performed at a first time point (forexample, prior to initiating therapy), and at a second time point (forexample, after one or more rounds of treatment), and a change a PKparameter between the first and the second time point provides a measureof the subject's response to treatment. Based on the results of themethods in this embodiment, an attending physician may recommend, forexample, (i) modifying the therapy to be administered to the subject(increased dosage, decreased dosage, increased frequency ofadministration, decreased frequency of administration, etc.), (ii)switching the subject to a different therapy, or (iii) maintaining acurrent course of therapy.

In other examples, treatment for ischemia, revascularization,inflammatory disease, infection, or other conditions may be evaluatedusing the method 1400. For instance, a change in a PK parameter maycorrespond to an increased or decreased blood perfusion in a tissue,which may be used to determine whether the tissue is responding to atreatment for ischemia, revascularization, inflammation, infection, orother condition. By way of non-limiting example, “inflammatory disease”refers to a disease or disorder characterized or caused by inflammation.“Inflammation” refers to a local response to cellular injury that ismarked by capillary dilatation, leukocytic infiltration, redness, heat,and pain that serves as a mechanism initiating the elimination ofnoxious agents and of damaged tissue. The site of inflammation mayinclude the lungs, the pleura, a tendon, a lymph node or gland, theuvula, the vagina, the brain, the spinal cord, nasal and pharyngealmucous membranes, a muscle, the skin, bone or bony tissue, a joint, theurinary bladder, the retina, the cervix of the uterus, the canthus, theintestinal tract, the vertebrae, the rectum, the anus, a bursa, afollicle, and the like. Such inflammatory diseases include, but are notlimited to, inflammatory bowel disease, rheumatoid diseases (e.g.,rheumatoid arthritis), other arthritic diseases (e.g., acute arthritis,acute gouty arthritis, bacterial arthritis, chronic inflammatoryarthritis, degenerative arthritis (osteoarthritis), infectiousarthritis, juvenile arthritis, mycotic arthritis, neuropathic arthritis,polyarthritis, proliferative arthritis, psoriatic arthritis, venerealarthritis, viral arthritis), fibrositis, pelvic inflammatory disease,acne, psoriasis, actinomycosis, dysentery, biliary cirrhosis, Lymedisease, heat rash, Stevens-Johnson syndrome, mumps, pemphigus vulgaris,and blastomycosis. Inflammatory bowel diseases are chronic inflammatorydiseases of the gastrointestinal tract which include, withoutlimitation, Crohn's disease, ulcerative colitis, and indeterminatecolitis. Rheumatoid arthritis is a chronic inflammatory diseaseprimarily of the joints, usually polyarticular, marked by inflammatorychanges in the synovial membranes and articular structures and by muscleatrophy and rarefaction of the bones.

The described systems and methods may also be used to aid in theprognosis of a subject suffering from cancer, ischemia,revascularization, inflammatory disease, infection, or other conditions.For example, the method 1500, alone or in any combination with one ormore of the blocks of method 100, may be used to determine a subject'sprognosis when blood perfusion is abnormal (e.g., more than normal inthe case of tumor angiogenesis, and less than normal in the case ofischemia, revascularization, or inflammation and infection). As usedherein, “prognosis” may, for example, refer to a likelihood that thesubject will or will not respond to a particular treatment regimen forthe relevant indication.

Initially, at block 1502, the method 1500 includes compare a determinedvalue of at least one PK parameter to a threshold or range. Thedetermined value may be a value determined using the method 100, forexample, or may be a value determined using other methods. In oneexample, the determined value may be a volume transfer constant, and thevolume transfer constant may be compared to a volume transfer constantthreshold or range. The comparison may determine, for instance, whetherthe determined value is within a predetermined range or a position ofthe determined value within the predetermined range (e.g., apercentile).

At block 1504, the method 1500 includes determine an effect of a therapyapplied to an interest region of tissue based on the comparison. In oneinstance, the method 1500 may be used to evaluate the early response ofa therapy, within days or hours of starting the therapy, for example.This may be advantageous relative to waiting for longer periods of timeto evaluate an effect, and then realizing that the therapy wasineffective.

In other examples, the described systems and methods may be used topredict the effect of a therapy applied to treat an abnormal conditionof blood perfusion, before the therapy is initiated. In still otherexamples, the described systems and methods may be used to evaluate theblood perfusion of patients who are at risk for some diseases, such ascancer. For example, individuals who have had a condition in the pastmay be at risk of developing a particular disease that can be evaluatedfor. As another example, a female patient with strong genetic markersfor breast cancer, who has dense breast tissue, smokes, and has ahistory of breast cancer in her family, may have enough risk factors tojustify screening the female patient for breast cancer.

The described applications are not meant to be limiting, and otherexamples are also possible. For instance, it is contemplated that theLRRM can be applied to the study of other phenomena. The derivation ofthe RAM is one application of LRRM. In line with this strategy, the LRRMcan also be applied to many other events which can also be characterizedby coupled ordinary differential equations, for example, quantitativemeasurement of enzymatic activity using hyperpolarized MRI, receptorbinding, and dynamic perfusion in non-biological systems.

In a further embodiment of all embodiments disclosed herein, the methodsmay further comprise treating the subject based on the results obtainedthrough carrying out the methods of the invention. As used herein.“treating” means accomplishing one or more of the following: (a)reducing or eliminating the relevant disorder, such as cancer, in thesubject; (b) reducing the severity of one or more symptoms of thedisorder; (c) limiting or preventing development of one or more symptomsof the disorder; (d) inhibiting worsening of one or more symptoms of thedisorder; and (e) limiting or preventing recurrence of one or moresymptoms of the disorder in subjects that were previously symptomaticfor the relevant symptom. The specific treatment to be provided to thesubject can be determined by an attending physician based on allfactors, including the relevant condition, the subject's specificcircumstances, and the results obtained using the methods of theinvention.

III. EXAMPLE SYSTEM

FIG. 16 is a simplified block diagram of an example system 1600 in whichthe present methods may be implemented. It should be understood thatthis and other arrangements described herein are set forth only asexamples. Those skilled in the art will appreciate that otherarrangements and elements (e.g., machines, interfaces, functions,orders, and groupings of functions, etc.) can be used instead and thatsome elements may be omitted altogether. Further, many of the elementsdescribed herein are functional entities that may be implemented asdiscrete or distributed components in conjunction with other componentsand in any suitable combination and location. Various functionsdescribed herein as being performed by one or more entities may becarried out by hardware, firmware, and/or software. For instance,various functions may be carried out by a processor executinginstructions stored in memory.

As shown in FIG. 16, example system 1600 may include a computing system1602 that is coupled to a display 1604. Generally, the computing system1602 may be configured to receive biomedical images 1606 or biomedicalimage data from another device (not shown). For instance, the computingsystem 1602 may receive the biomedical images 1606 from a database, aserver, another computing system, or a biomedical imaging device.Therefore, in some examples, the computing system 1602 may becommunicatively coupled to a biomedical imaging device that isconfigured to provide the biomedical images 1606. In other examples, thecomputing system 1602 may be configured to receive or access thebiomedical images from a database. Other examples may exist.

The biomedical images 1606 may take the form of any suitable type ofbiomedical imaging data. For example, the biomedical images 106 may bemagnetic resonance imaging (MRI) images such as dynamic contrastenhanced magnetic resonance imaging (DCE-MRI) images or ¹⁹F MRI images.In another example, the biomedical images may be optical images. Inother examples, the biomedical images may be ultrasonic images. Thus,the biomedical images may include any type of biomedical images orbiomedical image data.

The computing system 1602 may be configured to carry out the functionsdescribed herein. For example, the computing system 1602 may beconfigured to carry out any of the functions described with respect tothe example methods of FIGS. 1, 14, and 15 described above. As oneexample, the computing system 1602 may be configured to process thebiomedical images 1606 and determine one or more concentration-activitycurves based on the biomedical images 1606. For instance, an individualvoxel or individual pixel (or groups of voxels or pixels) may beanalyzed in one or more of the biomedical images 1606 to determine aconcentration of a contrast agent over time.

The contrast agent may be any type of substance used to enhance thecontrast of structures or fluids during medical imaging. For example,the contrast agent may be used to enhance the visibility of bloodvessels. For MRI examples, the contrast agent may be a substance thatalters the relaxation times of atoms within body tissues after oral orintravenous administration. As an example, the contrast agent may be aparamagnetic contrast agent, such as a gadolinium-based compound, thatalters the T₁ relaxation time of extracellular water. In anotherexample, the contrast agent may include fluorine. For instance, thecontrast agent may be a ¹⁹F nanoemulsion. For optical imaging examples,the contrast agent may include a colored dye. Other types of contrastagents are also contemplated.

The computing system 1602 may be further configured to fit theconcentration-activity curves to one or more pharmacokinetic models inorder to estimate at least one pharmacokinetic (PK) parameter. Forinstance, the computing system 1602 may be configured to determine aperfusion, permeability, and/or extracellular volume fraction by fittingtwo concentration activity-curves to a linear model based on applicationof a linear least square fitting (LLSQ) algorithm.

The computing system may also be configured to cause the display 1604 toprovide a visual indication of the value of at least one determined PKparameter. The display 1604 may be any type of graphical display deviceemploying one or more of any underlying display technologies (e.g.,cathode ray tube display, light-emitting diode display, liquid crystaldisplay, etc.).

The computing system 1602 of FIG. 16 may be implemented in, include, ortake the form of a computing device or components of a computing device,such as computing device 1702 shown in FIG. 17. As shown, computingdevice 1702 may include, without limitation, a communication interface1704, processor 1706, and data storage 1708, all of which may becommunicatively linked together by a system bus, network, and/or otherconnection mechanism 1714.

Communication interface 1704 typically functions to communicativelycouple computing device 1702 to other devices and/or entities. As such,communication interface 1704 may include a wired (e.g., Ethernet,without limitation) and/or wireless (e.g., CDMA and/or Wi-Fi, withoutlimitation) communication interface, for communicating with otherdevices and/or entities. Communication interface 1704 may also includemultiple interfaces, such as one through which computing device 1702sends communication, and one through which computing device 1702receives communication. Communication interface 1704 may be arranged tocommunicate according to one or more types of communication protocolsmentioned herein and/or any others now known or later developed.

Processor 1706 may include one or more general-purpose processors (suchas INTEL processors or the like) and/or one or more special-purposeprocessors (such as digital-signal processors or application-specificintegrated circuits). To the extent processor 1706 includes more thanone processor, such processors could work separately or in combination.Further, processor 1706 may be integrated in whole or in part withwireless-communication interface 1704 and/or with other components.

Data storage 1708, in turn, may include one or more volatile and/ornon-volatile storage components, such as magnetic, optical, or organicmemory components. As shown, data storage 1708 may include program data1710 and program logic 1712 executable by processor 1706 to carry outvarious functions described herein. Although these components aredescribed herein as separate data storage elements, the elements couldjust as well be physically integrated together or distributed in variousother ways. For example, program data 1710 may be maintained in datastorage 208 separate from program logic 1712, for easy updating andreference by program logic 1712.

Program data 1710 may include various data used by computing device 1702in operation. As an example, program data 1710 may include informationpertaining to biomedical image data and/or pharmacokinetic models.Similarly, program logic 1712 may include any additional program data,code, or instructions necessary to carry out the functions describedherein. For example, program logic 1712 may include instructionsexecutable by processor 1706 for causing computing device 1702 to carryout any of those functions described herein.

IV. EXAMPLE COMPUTER READABLE MEDIUM

In some embodiments, the disclosed methods may be implemented bycomputer program logic, or instructions, encoded on a non-transitorycomputer-readable storage media in a machine-readable format, or onother non-transitory media or articles of manufacture. FIG. 18 is aschematic illustrating a conceptual partial view of an example computerprogram product that includes a computer program for executing acomputer process on a computing device, arranged according to at leastsome embodiments presented herein.

In one embodiment, the example computer program product 1800 is providedusing a signal bearing medium 1802. The signal bearing medium 1802 mayinclude one or more programming instructions 1804 that, when executed byone or more processors may provide functionality or portions of thefunctionality described herein. In some examples, the signal bearingmedium 1802 may encompass a computer-readable medium 1806, such as, butnot limited to, a hard disk drive, a Compact Disc (CD), a Digital VideoDisk (DVD), a digital tape, memory, etc. In some implementations, thesignal bearing medium 1802 may encompass a computer recordable medium1808, such as, but not limited to, memory, read/write (R/W) CDs, R/WDVDs, etc. In some implementations, the signal bearing medium 1802 mayencompass a communications medium 1810, such as, but not limited to, adigital and/or an analog communication medium (e.g., a fiber opticcable, a waveguide, a wired communications link, a wirelesscommunication link, etc.). Thus, for example, the signal bearing medium1802 may be conveyed by a wireless form of the communications medium1810. It should be understood, however, that computer-readable medium1806, computer recordable medium 1808, and communications medium 1810 ascontemplated herein are distinct mediums and that, in any event,computer-readable medium 1808 is a physical, non-transitory,computer-readable medium.

The one or more programming instructions 1804 may be, for example,computer executable and/or logic implemented instructions. In someexamples, a computing device such as that shown in FIG. 17 may beconfigured to provide various operations, functions, or actions inresponse to the programming instructions 1804 conveyed to the computingdevice by one or more of the computer readable medium 1806, the computerrecordable medium 1808, and/or the communications medium 1810.

The non-transitory computer readable medium could also be distributedamong multiple data storage elements, which could be remotely locatedfrom each other. The computing device that executes some or all of thestored instructions could be a computing device, such as the computingdevice illustrated in FIG. 17. Alternatively, the computing device thatexecutes some or all of the stored instructions could be anothercomputing device.

V. CONCLUSION

It is intended that the foregoing detailed description be regarded asillustrative rather than limiting and that it is understood that thefollowing claims including all equivalents are intended to define thescope of the invention. The claims should not be read as limited to thedescribed order or elements unless stated to that effect. Therefore, allembodiments that come within the scope and spirit of the followingclaims and equivalents thereto are claimed as the invention.

We claim:
 1. A computer-implemented method, comprising: determining,based on a set of biomedical images, a first concentration-activitycurve and a second concentration-activity curve, wherein each of thefirst concentration-activity curve and the second concentration-activitycurve indicates a concentration of a respective contrast agent within arespective region of tissue; determining a value of at least onepharmacokinetic (PK) parameter based on the first concentration-activitycurve and the second concentration-activity curve, wherein determiningthe value of the at least one PK parameter comprises: (i) determining alinear model that relates the first concentration-activity curve to thesecond concentration-activity curve, wherein the linear model comprisesthe at least one PK parameter; and (ii) determining the value of the atleast one PK parameter based on application of a linear least squarefitting (LLSQ) algorithm to the linear model; and causing a graphicaldisplay to provide a visual indication of the value of the at least onePK parameter.
 2. The computer-implemented method of claim 1: wherein thefirst concentration-activity curve indicates a concentration of a firstcontrast agent in an interest region of tissue; and wherein the secondconcentration-activity curve indicates a concentration of the firstcontrast agent in a reference region of tissue.
 3. Thecomputer-implemented method of claim 1, wherein the first contrast agentcomprises a paramagnetic contrast agent.
 4. The computer-implementedmethod of claim 1: wherein the first concentration-activity curveindicates a concentration of a first contrast agent in an interestregion of tissue; and wherein the second concentration-activity curveindicates a concentration of a second contrast agent in the interestregion of tissue.
 5. The computer-implemented method of claim 4, whereineach of the first contrast agent and the second contrast agent comprisefluorine.
 6. The computer-implemented method of claim 4: wherein thefirst contrast agent comprises a first color; and wherein the secondcontrast agent comprises a second color.
 7. The computer-implementedmethod of claim 4: wherein the first contrast agent is responsive to abiomarker; and wherein the second contrast agent is unresponsive to thebiomarker.
 8. The computer-implemented method of claim 1, wherein the atleast one PK parameter comprises one or more of the following: a volumetransfer constant (K^(trans)) and an extravascular extracellular volumefraction (v_(e)).
 9. The computer-implemented method of claim 1, whereinthe set of biomedical images comprises a set of magnetic resonanceimages.
 10. The computer-implemented method of claim 1, wherein the setof biomedical images comprises a set of optical images.
 11. Thecomputer-implemented method of claim 1, further comprising determining aresponse-indication corresponding to a treatment of an interest regionof tissue based on the value of the at least one PK parameter.
 12. Thecomputer-implemented method of claim 11, wherein the treatment comprisesone or more of the following: an anti-angiogenic therapy, a treatmentfor ischemia, a treatment for revascularization, and a treatment forinflammation and infection.
 13. The computer-implemented method of claim1, further comprising determining that an interest region of tissue is amalignant tumor based on the value of the at least one PK parameter. 14.A non-transitory computer-readable medium having instructions storedthereon, the instructions comprising: determining, based on a set ofbiomedical images, a first concentration-activity curve and a secondconcentration-activity curve, wherein each of the firstconcentration-activity curve and the second concentration-activity curveindicates a concentration of a respective contrast agent within arespective region of tissue; determining a value of at least onepharmacokinetic (PK) parameter based on the first concentration-activitycurve and the second concentration-activity curve, wherein determiningthe value of the at least one PK parameter comprises: (i) determining alinear model that relates the first concentration-activity curve to thesecond concentration-activity curve, wherein the linear model comprisesthe at least one PK parameter; and (ii) determining the value of the atleast one PK parameter based on application of a linear least squarefitting (LLSQ) algorithm to the linear model; and causing a graphicaldisplay to provide a visual indication of the value of the at least onePK parameter.
 15. The non-transitory computer-readable medium of claim14: wherein the first concentration-activity curve indicates aconcentration of a first contrast agent in an interest region of tissue;and wherein the second concentration-activity curve indicates aconcentration of the first contrast agent in a reference region oftissue.
 16. The non-transitory computer-readable medium of claim 14:wherein the first concentration-activity curve indicates a concentrationof a first contrast agent in an interest region of tissue; and whereinthe second concentration-activity curve indicates a concentration of asecond contrast agent in the interest region of tissue.
 17. A systemcomprising: at least one processor; a computer-readable medium; andprogram instructions stored on the computer-readable medium andexecutable by the at least one processor to: determine, based on a setof biomedical images, a first concentration-activity curve and a secondconcentration-activity curve, wherein each of the firstconcentration-activity curve and the second concentration-activity curveindicates a concentration of a respective contrast agent within arespective region of tissue; and determine a value of at least onepharmacokinetic (PK) parameter based on the first concentration-activitycurve and the second concentration-activity curve, wherein determiningthe value of the at least one PK parameter comprises: (i) determining alinear model that relates the first concentration-activity curve to thesecond concentration-activity curve, wherein the linear model comprisesthe at least one PK parameter; and (ii) determining the value of the atleast one PK parameter based on application of a linear least squarefitting (LLSQ) algorithm to the linear model.
 18. The system of claim17: wherein the first concentration-activity curve indicates aconcentration of a first contrast agent in an interest region of tissue;and wherein the second concentration-activity curve indicates aconcentration of the first contrast agent in a reference region oftissue.
 19. The system of claim 17: wherein the firstconcentration-activity curve indicates a concentration of a firstcontrast agent in an interest region of tissue; and wherein the secondconcentration-activity curve indicates a concentration of a secondcontrast agent in the interest region of tissue.
 20. The system of anyof claim 17, further comprising a dynamic contrast enhanced magneticresonance imaging (DCE-MRI) device communicatively coupled to thesystem, wherein the DCE-MRI device is configured to provide the set ofbiomedical images to the system.